Astronomy & Astrophysics manuscript no. 
(will be inserted by hand later) 



Nonlinear effects in time-resolved spectra of DAVs 

J. Ising and D. Koester 

Institut fur Theoretische Physik und Astrophysik, Universitat Kiel, D-24098 Kiel, Germany 
(ising,koester@astrophysik.uni-kiel.de) 

February 1, 2008 

Abstract. Numerical simulations of light curves of variable DA white dwarfs (ZZ Ceti stars) predict flux amplitudes 
with surface distributions different from the spherical harmonics of the pulsation mode in deeper layers. In contrast 
to the results of the perturbation analysis by Goldreich and Wu this is also true for the fundamental period of 
the flux variation. As a consequence normalized amplitude spectra depend not only on the mode number I but 
also on pulsation amplitude and inclination. Another new result is that with increasing amplitude of the pressure 
variation below the convection zone the flux variation at the surface goes through a maximum and then decreases 
again. 
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1. Introduction 

The ZZ Ceti stars are the coolest class of pulsating white 
dwarfs. They are multi-periodic g-mode pulsators with pe- 
riods in the range of 100s to 1000s. In contrast to the 
PG 1159 objects and the variable DB only a few modes 
are observed in each individual ZZ Ceti star. The num- 
ber of observed pulsation periods, the amplitudes and the 
shape of the light curves of these stars change syst emat i- 
cally within the instability strip (see e.g. Clemens |l993|) . 
At the hot (blue) edge small amplitude pulsations are ob- 
served with sinusoidal light curves. The cooler ZZ Ceti 
stars show more modes with large amplitudes of more than 
10% with overtones and additional frequencies, formed by 
sums and differences of the fundamental frequencies of the 
modes. 

The small number of unstable modes forms an obsta- 
cle to standard mode identification techniques (identifying 
the "quantum numbers" I and m of the spherical harmon- 
ics of a mode), which rely mostly on periods and period 
spacings of unstable frequencies in the power spectrum. In 
recent years time-resolved spectroscopy of white dwarfs 
has become possible and was used to identify the spheri- 
cal harmonics of the modes. This mode identification tech- 
nique uses the wavelength dependence of normalized pul- 
sation amplitudes and takes advantage of the wavelength 
dependence of the limb darkening and different cancella- 
tion of the flux variation for different spherical harmonics 
(Robinson et al. 1982; Brassard et al. 1995; Robinson et 
al. |1995|; Kepler et al. 2000| ; Clemens and Kerkwijk p000| ) . 



This diagnostic technique depends on the assumption that 
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the flux amplitude on the visible surface varies with the 
spherical harmonic (I, m) of the mode. 

This assumption — which we will call "linear" - is 
subject to criticism, at least for larger amplitudes. For 
large amplitudes non-sinusoidal light curves are observed, 
in contrast to the sinusoidal light curves for small ampli- 
tudes. Because the ZZ Ceti stars are non-radial pulsators, 
the amplitude varies over the surface of the star and one 
would expect a variation of the shape of the light curves 
over the surface for large amplitudes as well. This is in 
contradiction to the assumption of a flux variation with 
spherical harmonics. For a quantitative description of the 
integrated light curve and spectra we will need a theo- 
retical understanding of the reaction of the local surface 
flux to an imposed pressure variation at deeper layers. 
This theory should demonstrate how the surface variation 
becomes increasingly non-sinusoidal with increasing am- 
plitude of the pressure variation. 

Such a description was provided by Brickhill (1983), 
who proposed, that the thin convection zone below the 
surface of the ZZ Ceti stars has an important effect on 
the light curves. Brickhill showed in a series of papers, that 
numerical simulations of the light curves can qualitatively 
reproduce the observed flux variations. This result was 
confirmed by Wu ( |1998|) and Wu ( [20001 ). 

Gol dreich & Wu ( |l999| , pJ99a| ) and Wu & Goldreich 
(1998a) studied the driving mechanism of the ZZ Ceti 
stars in a linear perturbation analysis. Their work lead to 
an impressive progress in the understanding of the driv- 
ing of the pulsation modes, but it their linear analysis they 
consider only the entr opy change in a convection zone of 
constant depth. Wu (1998,2000) extended this work to 



2 



J. Ising and D. Koester: Nonlinear effects in time-resolved spectra of DAVs 



take into account the depth changes, but used only static 
models of convection zone. 

The numerical integrations we have performed for this 
paper show that for larger amplitudes the depth varies 
significantly during one cycle. At each time step during the 
pulsation the structure remains different from any static 
structure, leading to modifications of the nonlinear effects, 
which appear already in the perturbation analysis using 
static structures, as described in this paper. 

We extend Brickhill's work and show that the con- 
vection zone strongly influences the spatial distribution 
of surface amplitudes, which for larger amplitudes cannot 
anymore be described by spherical harmonics. This can- 
not be neglected in the technique of mode identification 
by time-resolved spectroscopy. In the next section we de- 
scribe the model used for our numerical calculations. First 
we will introduce the relatively simple system of ordinary 
differential equations to solve the time-dependent struc- 
ture of the envelope of ZZ Ceti stars, obtained by allowing 
only a restricted time dependence of the basic equations. 
To examine the validity of this restriction for large am- 
plitude pulsation we introduce the full time-dependent 
equations in the following subsection. The boundary con- 
ditions have a large influence on the results of the numeri- 
cal simulations; the following subsection is devoted to this 
topic. The 3. section describes the results of a light curve 
simulation for a plane-parallel column with a fixed pul- 
sation amplitude. In the 4. section we consider that the 
pulsation amplitude varies over the surface for non-radial 
pulsators like the ZZ Ceti stars, and describe the surface 
distribution of the flux for different maximum amplitudes 
and spherical harmonics. The technique used to calculate 
the total flux together with the expected consequences for 
the time-dependent spectroscopy are discussed in the 5. 
section. Finally we give our main conclusion in section 6. 

2. The numerical model 



Within her analysis Wu ( 1998 200C| calculates the local 
response of the surface flux to a sinusoidal variation of 
pressure and flux in a layer at the bottom of the convection 
zone. Due to the change in extent and structure of the 
convection zone its reaction to the perturbation of the 
input flux is nonlinear, leading to the appearance of higher 
harmonics in the surface flux (Wu 2000)). 

The main result is that the amplitude of the funda- 
mental period depends linearly and that of the first over- 
tone quadratically on the amplitude below the convection 
zone. The phase delays between surface flux and pressure 
variation are constant, independent of amplitude. This 
important result guarantees that the spatial distribution 
over the stellar surface for the fundamental period is still 
given by the same spherical harmonic, which describes the 
variation in the deeper layers. This justifies the usual as- 
sumptions made in the analysis of time-dependent spec- 
troscopy. 

The starting point of our own study is the numerical 
model of Brickhill. In this model the convective layer is 



considered as a plane parallel column. The lower bound- 
ary of this column lies below the surface convection zone, 
the upper boundary is directly below the photosphere. The 
sinusoidal relative pressure variation is assumed to be con- 
stant through the whole column, an assumption justified 
in the work by Brickhill and Goldreich and Wu. 

In the simulation we calculate the transfer of energy 
and the time-dependent temperature structure numeri- 
cally in a way similar to Brickhill's work. Some improve- 
ments over his work arc state-of-the-art equation of state 
from Saumon ct al.( 1995) and OPAL opacity data (Iglesias 
and Rogers, 1996), and our own model atmosphere grids, 
which are described in Finley et al. (1997). 

In this paper we want to show the principal effects for 
the analysis of time-resolved spectroscopy. Since we are 
not aiming at a detailed comparison with observations, 
it is sufficient to investigate only one equilibrium model 
with one effective temperature and one pulsation period. 
We take a stellar model with T e g = 11350K and log<? = 8 
and a pulsation period of 100s. We use only one parameter 
for the mixing length theory (ML2/a = 0.6); the thermal 
relaxation time scale is 4.5 s. With this choice our model 
is representative for the conditions near the blue edge of 
the instability strip. All our calculations and results of this 
study are based on this set of parameters. 

2.1. The basic equations 

The main aim of the model calculation is to determine the 
temperature structure as a function of time. The tempe- 
rature disturbance ST for a static model can be calculated 
from 



ST = 



SQ 



C P AM 



V ad T- 



.SP 
P~' 



(1) 



Equation (Q) is equivalent to the first law of thermody- 
namics (e.g. Cox 198C, p. 37) and can be found in Brickhill 
( 1983[ ) as well. The first term on the right describes the 
non-adiabatic heating by the additional heat energy SQ of 
a mass AM with the heat capacity C p . The second term is 
the adiabatic contribution to the temperature change due 
to the relative pressure variation ^p- in an environment 
with an adiabatic temperature gradient V a d- 

To calculate the heat energy SQ we consider a finite 
volume element with the cross section A and a height Az. 
Introducing a discretization and finite volume elements 
at this point, we can avoid the use of partial differential 
equations. In a plane-parallel symmetry the flux has only 
a vertical component. If AF is the difference between the 
outgoing and the incoming flux for a mass element, the 
heat energy accumulated in a time interval At is given by 



At 

dtAAF 

or written in differential form 
dSQ 



dt 



= AAF 



(2) 



(3) 
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To calculate the mass AM of the volume element we 
follow the procedure of Brickhill. We take the pressure of 
the non-pulsating star in hydrostatic equilibrium to define 
the vertical depth scale. Following Brickhill we introduce 
Q = In P as the independent variable of the problem. The 
mass element for a fixed pressure difference Ad> can be 
written as 



AM = A 



A 



dz g(t) 



U>1 



Gj\ -\-ACj 



A 



wi+Au P(t) 

duj . 

a>i 9 



(4) 



As B rickh ill assumed and was confirmed by Goldreich & 
Wu ( 1999 ) the relative pressure variation can be taken as 
constant through the whole column. We further assume 
a constant surface gravity g throughout the column. This 
simplifies the expression for the mass element: 



AM = A 



exp(f (t)) 



exp(w) / dtij' exp(o)') 



with 



f(x) = cxp(x) - 1. 

In a first order approach we can combine equation (|l|) , 
(||) and (0), by taking the derivative with respect to the 
time of (fy). In this first order approach we are only con- 
cerned with the time-dependency of ST, SQ and ^p-. We 
get 



dT _ AFg 
~dt ~ G p Pf(AQ) 



d 5P 



(6) 



For all the mass elements of the vertical column, equa- 
tion (^) defines a system of ordinary differential equations 
(ODE) for the temperature structure. The fluxes describe 
the coupling of the equations at each grid-point of the col- 
umn. In our approximation we calculate the flux locally 
from the actual temperatures and temperature gradients. 
It is given by 



F = F r + F c 



(7) 



where the radiative flux F r can be calculated with the dif- 
fusion approximation and the convective flux F c by means 
of the mixing length theory. We use a time-independent 
formulation of the mixing length theory and assume an 
instantaneous adjustment of the total flux to the temper- 
ature and pressure structure of the column. To check the 
validity of this assumption we have also used the time- 



dependent Kuhfuss model for convection (Kuhfuss, 1986; 
Wuchterl & Feuchtinger, 1998) and found no differences 



between the time-dependent and the time-independent 
model calculations. This result holds strictly, if the ki- 
netic energy of the convection does not contribute to the 
thermal energy, but the differences are also small, if we 
take such a coupling into account. 



The rapid variation of the flux with the temperature 
leads to short relaxation time scales near the upper bound- 
ary of the column. This leads to a stiff behavior of the 
ODE. Consequently we use a Gear solver (Gear, 1971) to 
integrate equation (^). 

2.2. The full time-dependent equations 

The differential formulation of (|l|) leads in principle to 
reliable solutions for large perturbations as well. But to 
obtain equation ([|) we have ignored all time-dependencies 
except for ST, SQ and . For large amplitudes {^S- > 5%) 
this cannot be justified any longer, since temperature and 
specific heat capacity can change by large amounts. In a 
first step to correct for this we can replace all quantities 
in equation (|^) by their current values instead of those 
of the equilibrium model. Since this may not be sufficient 
for short pulsation periods (~ 100s), we have modified 
equation (|^) and taken into account all terms produced 
by the derivative of (m with respect to time 



dT 



(5) [l + x 1 -x 2 )- = 



AFg T d_5P_ 

C p Pf(A*) +Vadl dt P 



■x 3 + x 4 (8) 



In the following all thermodynamic quantities are func- 
tions of T and P, and a partial derivative with respect to 
T means that P is to be kept constant and vice versa. The 
correction terms x\ and X3 due to the non-adiabatic term 
of equation (m are 



Xi 



X3 



SQ dC p 
CI AM dT 

SQ dC p dP 
CfAM~dP~dt 



The correction terms X2 and £4 are produced by the 
derivative of the adiabatic term of (0) 



X2 



X4 



a< 1 



T 



' ad 



V dT 
gVad T SP dP 
dP P dt 



5P 
~P 



The derivative of the mass element with respect to 
the time is calculated again under the assumption, that 
the relative pressure perturbation is constant through the 
whole column. Then this variation is inversely related to 
the change of the horizontal area 

5 A _ SP 
~A ~ ~~P ' 
This leads to 



dAM fdA^, x Arx dP 



/(A£) 
9 



0. 



(9) 



In equation (]8j) appear both the time-derivative of SQ 
and SQ itself. This changes the structure of the equations 
from ordinary differential equations to a set of integro- 
differential equations. We solve this system by a splitting 
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procedure. In the inner step we solve the differential equa- 
tions for a time step short compared to the pulsation pe- 
riod with constant SQ. After this time step SQ is updated. 
For time steps of about 0.25% of the period the achieved 
accuracy is high enough, since the integral only appears 
in higher oder terms. 

The contributions of x± and xi are important for 
large pulsation amplitudes, caused by the large and rapid 
changing of the temperature, whereas all terms including 
a derivative of V a d are unimportant. Table ([!]) compares 
the amplitudes of the fundamental and the first overtone 
of the photospheric flux variation as a result of (||) with 
(PD for our model with a period of 100s. For (g) we take 
the current values of all quantities instead of the values of 
the equilibrium model. 



For the lower boundary we assume the adiabatic ap- 
proximation. The matter in the region below the convec- 
tion zone is well approximated as a completely ionized 
hydrogen plasma. Under this assumption the flux pertur- 
bation is given by (see Goldreich & Wu, p99|) 



SF _ 3 - 3k b - 2k t SP 
T ~ 5 ~P~ ' 



(10) 



where n g and kt are the logarithmic derivatives of the 
opacity with respect to the density and temperature. This 
adiabatic assumption is the weakest point of the approx- 
imations in the numerical model. In order to test this as- 
sumption and also to find the optimum location of the 
lower boundary we have made a numerical experiment. 





ODE 


IDE 




Ax[%] 


A 2 [%] 


A x [%] 


A 2 [%] 


0.02 


1.71 


0.24 


1.71 


0.25 


0.04 


3.64 


1.17 


3.65 


1.23 


0.06 


5.14 


2.01 


5.18 


2.13 


0.08 


6.17 


2.61 


6.29 


2.83 


0.10 


6.78 


2.95 


7.04 


3.32 


0.12 


6.97 


2.97 


7.45 


3.56 


0.14 


6.74 


2.67 


7.54 


3.54 


0.16 


6.09 


2.01 


7.31 


3.25 


0.18 


5.17 


1.17 


6.78 


2.68 



Table 1. Fourier amplitudes for the photospheric flux 
variation. Column 2 and 4 display the amplitude of the 
fundamental for equation (||) (ODE) and equation (|J) 
(IDE), respectively. Column 3 and 5 show the amplitudes 
of the first overtone. All quantities are functions of the 
pressure amplitude ^p-. 
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Fig. 1. The variation of the photospheric flux (solid line) 
and the flux variation at the lower boundary (dashed line) 
are plotted for a 1% (upper panel) and a 5% pressure 
amplitude (lower panel). 



The differences for large pressure amplitudes are signif- 
icant in the solutions of (^) and (||), but the qualitative be- 
havior is the same in both cases. Since the numerical effort 
is not very much larger for solving the IDE this method 
will be implemented in future versions of our code. The 
ODE solution (||) is adequate for a qualitative discussion 
and all results given in this paper are based on (|^) with 
current values for all quantities. One should remember, 
however, that the results for large amplitudes — though 
qualitatively correct — may be affected by this approxi- 
mation. 

2.3. Boundary conditions 

The basic system of ODE (||) looks like an initial value 
problem, but this is only partially correct. To obtain equa- 
tion (^|) we have introduced a discretization of the column, 
to replace the spatial derivatives of the flux with finite dif- 
ferences. The mathematical structure of a partial differen- 
tial equation survives in the coupling of the ODEs and the 
necessity to define boundary conditions for the flux. 



We assume a column with a lower boundary much 
deeper than the bottom of the convection zone with a con- 
stant pressure amplitude in the whole column. The depth 
of the lower boundary is limited by the growing comput- 
ing costs, because a pulsation simulation has to be done 
for at least one thermal time scale of the whole column to 
reach a stationary situation. For our model we choose a 
thermal time scale of 1500s and a value of ui = 22 to fix the 
lower boundary for this experiment. At this lower bound- 
ary we apply a sinusoidal flux perturbation and study the 
resulting flux at a point clearly below the maximum ex- 
tent of the convection zone, but far away from the lower 
boundary. At this point the flux is no longer sinusoidal 
but shows a fundamental with amplitude Ai and a first 
overtone with amplitude A2. The phases are shifted rela- 
tive to the bottom flux by 4>\ and 4>2 respectively, where a 
positive phase means that the variation is lagging behind 
the perturbation. One should note that these effects would 
be absent in a completely radiative star; they are caused 
by the delayed reaction of the convection zone, which de- 
termines the outer boundary for the envelope structure. 
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Tabic g shows the results at the point with Co — 19.45. 
The first overtone at this point remains at least one or- 
der of magnitude smaller than the fundamental amplitude 
and the phase shift between pressure and flux variation re- 
mains small. We follow Brickhill in our conclusion that the 
adiabatic flux variation at the lower boundary is a good 
approximation if the propagation zone of the g-modes is 
far away from the convective layer. This is true for short 
period pulsations as considered in our model, which rep- 
resents the "blue edge" of the instability strip. For our 
pulsation simulations we then fixed our lower boundary 
at Co = 20. 



p 




A x [%] 


0i n 


A 2 [%] 


fa [°] 


0.02 


4.40 


4.72 


2.67 


0.05 


-98.43 


0.04 


8.80 


9.43 


2.67 


0.19 


-97.35 


0.06 


13.20 


14.09 


2.84 


0.40 


-97.80 


0.08 
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3.12 


0.69 


-97.80 


0.10 
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23.05 


3.43 


1.08 


-97.52 


0.12 


26.40 


27.39 


3.75 


1.56 


-97.35 


0.14 


30.80 


31.59 


4.09 


2.10 


-97.92 


0.16 


35.20 


35.61 


4.47 


2.72 


-99.01 


0.18 


39.60 


39.45 


4.88 


3.44 


-100.32 



Table 2. Influence of non-adiabatic terms at the lower 
boundary: Results of the experiment described in the text. 
The 2 nd column displays the amplitude of the sinusoidal 
flux variation at the lower boundary. The 3 rd and 4 th col- 
umn is the amplitude and phase of the fundamental at 
the point Co = 19.45, which is always below the convection 
zone. The first overtones are shown in column 5 and 6. 



To find an upper boundary condition we make the as- 
sumption, that the photosphere reaches the static struc- 
ture for each time step instantaneously. As discussed e.g. 
by Brickhill (1983) this approximation is valid for small 
thermal time scales compared with the pulsation period. 
We choose the upper boundary at the pressure point with 
the optical depth « 10, where this condition is clearly 
fulfilled. The advantage to take the optical depth clearly 
greater than 1 at the upper boundary of the column is, 
that the time-dependent calculations can be done with the 
diffusion approximation for the radiative flux. The deter- 
mination of the photospheric structure requires a detailed 
radiative transfer calculation, but without consideration 
of terms with time-derivatives. To connect these two dif- 
ferent calculations we assume, that the combination of the 
actual flux, pressure and temperature is the same at the 
upper boundary of the column as in one specific static pho- 
tospheric model. To find the flux we use a grid of static 
models and find the known combination of pressure and 
temperature in this grid for each time step separately and 
interpolate the flux from these models. This procedure 
avoids a further linearization at the upper boundary, but 
is in principle equivalent to Brickhill's method. 



100 




500 



Fig. 2. Shape of the flux variation for different pressure 
amplitudes. The bottom line is the light curve for 1% am- 
plitude, the increment between curves is 1%. The light 
curve at the top of the panel is for 20% pressure ampli- 
tude. 



3. Numerical results for a column 

In this section we study the reaction of the convection zone 
and the shape and amplitude of the surface flux variation 
for a column on the stellar surface. 

For small amplitude pulsations with 1% pressure vari- 
ation (see upper panel in Fig. |l]) the photospheric flux 
(solid line) is shifted in phase relative to the flux of the 
adiabatic layers (dashed line) and the amplitude is re- 
duced. As shown by Brickhill ( 1983| ), the phase shift and 
the reduced amplitude is caused by the delay in the heat 
transfer due to the necessary adjustment of the convective 
layers to a change in the input flux. If the flux is increased, 



6 



J. Ising and D. Koester: Nonlinear effects in time-resolved spectra of DAVs 



0.08 


p.a/f 


_ fundamental 


0.06 


jy I /p-ct. 




0.04 


_ / p. a./ / 






. / // 1 ! 


overtone 


0.02 


/ .... ■ 2 nd 


overtone-.. ^ . 


0.00 


// ' ' 





0.00 



0.05 



0.10 
dP/P 



0.15 



1.00 



A 0.98 



0.96 



0.94 



0.00 




Fig. 3. Amplitudes of the Fourier coefficients relative to 
the mean flux for a local column. The numerical results 
are labeled with the name of the harmonic. The curves 
marked with p. a. indicate the linear and quadratic rela- 
tions predicted by the perturbation analysis. 
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Fig. 4. Phases of the fundamental and overtones relative 
to the pressure variation in a local column for the nu- 
merical results. The linear perturbation analysis predicts 
constant phases with values, which depend on the equilib- 
rium model. 



the depth of the convective layer is shrinking as required 
by a static envelope solution corresponding to a higher ef- 
fective temperature. A large amount of heat is necessary to 
change convective into radiative layers and to increase the 
heat content of the remaining convection zone, which has 
to be supplied by the input flux. Consequently, the pho- 
tospheric flux follows the input flux with a delay of the 
order of the thermal time scale of the convective layer. If 
the input flux varies periodically with a period compara- 
ble to the thermal time scale, this delay results in a phase 
shift and an amplitude reduction. 

Fig. ^ summarizes this variation of the light curves 
with changing pressure amplitudes. 



Fig. 5. Time-averaged flux relative to the flux below the 
convection zone. 

For small amplitudes the variation of the surface flux 
remains sinusoidal and the amplitude increases linearly 
with the pressure amplitude as predicted by the pertur- 
bation analysis. 

For larger amplitudes, e.g. 5% (lower panel Fig. 0) the 
extent of the the convective zone varies significantly dur- 
ing the pulsation cycle, and the stratification may even 
become completely radiative shortly after maximum com- 
pression. At that instant the photospheric flux follows di- 
rectly the adiabatic flux variation, causing the steep in- 
crease and sharp maxima in the light curves. 

Increasing the pressure amplitudes even further (20%), 
the surface flux variation becomes sinusoidal again, with 
small amplitudes. In this CclS6, cL large fraction of the heat 
flux is transformed into mechanical energy, and the re- 
maining average flux corresponds to a model with lower 
(time-averaged) effective temperature and thicker convec- 
tion zone. The thermal time constant of this convection 
zone becomes large compared to the pulsation period, re- 
sulting in a strong reduction of the flux amplitude. 

We emphasize that this conversion of heat to me- 
chanical energy in the column (which remains qualita- 
tively the same with the higher numerical accuracy of 
the IDE method) does not directly tell us anything about 
the global excitation or damping of the mode. Our calcu- 
lation is based on energy arguments in an open system: 
the enforced variations of the pressure and column cross 
section carry away any difference of incoming and outgo- 
ing heat flux. In order to really study the excitation of 
the mode one would need to calculate the behavior of the 
mode throughout the whole star. It is quite possible that 
the increase of the mechanical energy in the outer enve- 
lope is more than balanced by a damping in the interior; 
it is also possible that the energy goes into the excitation 
of other modes. This problem will be studied further in 
future simulations, which include the study of the P dV 
term in Eq. [1} 

The sharp maxima in the light curves imply the ap- 
pearance of higher harmonics in the Fourier decomposi- 
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tion. Fig. |3| demonstrates, how the relative amplitudes of 
fundamental, first, and second overtone increase with in- 
creasing pressure amplitude. For low amplitudes the re- 
sults confirm the predictions of the perturbation analy- 
sis (linear and quadratic dependence for fundamental and 
first overtone); for amplitudes larger than 5%, however, 
strong deviations become apparent. 

The sign and even the magnitude of the phase shifts in 
Fig. |] approximately agree with the predictions of the per- 
turbation analysis: the fundamental for the flux is lagging 
behind the pressure, whereas the in the first harmonic it 
is leading. However, quite unexpectedly from the predic- 
tions of Goldreich and Wu, for very small amplitudes the 
phases do not become constant. The shift demonstrates 
the change of the relative contributions from the com- 
pression in the convection zone, which is strictly in phase, 
and the flux variation below the convection zone, which 
is phase shifted during the propagation through the con- 
vective layer. While at very low amplitudes the numerical 
result may be questionable, we are confident that this is 
not the case for pressure amplitudes larger than about 
1%. This result needs further study, and we will attempt 
in a future paper to reconcile the results of numerical and 
perturbational analysis. 

Another difference to the perturbation analysis of 
Goldreich and Wu is the dependence of the mean local 
flux on the pressure amplitude (see Fig. ^). The flux from 
below the convection zone is reduced by a fraction of en- 
ergy that is converted into mechanical energy to drive the 
pulsation mode. This fraction grows quadratically with 
the pressure amplitude. Consequently the mean flux is re- 
duced very efficiently for large amplitudes. 

4. Surface distribution 

In the preceding sections we have discussed the local reac- 
tion of the photosphcric flux to an assumed fixed pressure 
amplitude. In this section we consider the whole surface of 
the star, assuming that the spatial structure of the pres- 
sure perturbation as well as the flux perturbation below 
the convection zone is described by a spherical harmonic 
of mode numbers I and m. 

For small amplitudes — because of the linear depen- 
dence of surface amplitudes on pressure amplitudes — the 
distribution of the surface flux amplitude follows directly 
the pressure and is therefore described by the same spheri- 
cal harmonic. But this is only true if one ignores the small 
phase changes of the flux variation. Due to this phase shift 
the maximum flux amplitude is reached at a later time for 
smaller amplitudes. This effect leads to a deviation of the 
photospheric flux distribution from the spherical harmonic 
of the mode, but the flux amplitude at the surface is still 
a monotonous function of the pressure amplitude. This 
changes with increasing amplitude: because of the non- 
monotonous local relation the maximum flux amplitude 
at the surface does not coincide with the maximum of the 
pressure amplitude and e.g. for m = modes the maxi- 
mum flux amplitude is not reached at the poles. Since also 
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Fig. 6. Surface distribution of the effective temperature 
near the maximum of the integrated flux. The gray-level 
of the plot (top, for an I = 2 mode) is normalized to 
the maximum T e g (white) and minimum T e ft (black) on the 
surface. The middle panel shows in graphical form the 
temperature distribution with angle /i = cos $ for I = 1 
at three different times, the bottom panel the same but 
for 1 = 2. The maximum pressure amplitude is 2% (small 
amplitude) . 



the phases change with amplitude, different latitudes on 
the surface reach the maximum flux at different times. 

Fig. U to H show the surface distribution for I = 
1,2, m = modes as three-dimensional plot with gray- 
level indicating the effective temperature, and in graphical 
form the dependence on polar angle \i = cosi9. Shown are 
three different cases for small, medium, and large pressure 
amplitudes. 
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Fig. 7. Similar to Fig. |5|, but for a maximum pressure 
amplitude of 5% 
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Fig. 8. Similar to Fig. |[ but for a maximum pressure 
amplitude of 20% 



For I = 2 modes the change due to the nonlinear rela- 
tion between pressure amplitude and surface flux ampli- 
tude may be quite dramatic: whereas for small amplitudes 
the maximum amplitudes are at the poles, with a smaller 
relative maximum at the equator, the situation can be 
completely reversed at large amplitudes (see Fig. |8|). 

5. Flux integration and time-resolved spectra 

From these results it follows that the surface distribution 
of the flux variation in general has a much more compli- 
cated form than is assumed traditionally (meaning spher- 
ical harmonics describing all variations on the stellar sur- 
face). In this section we want to answer two questions: 
First we want to find an optimal strategy to calculate the 
wavelength dependent amplitude spectra, and second we 
want to discuss the major effects on time-resolved spec- 



tra one has to expect, if the linear assumption is not any 
longer appropriate. 

The general problem is to find the time-dependent flux 
T\ = T\{yi,t), where ft is the direction to the observer. 
To integrate the flux we need the specific intensity at each 
point of the surface at each time in the direction ft. In gen- 
eral, the specific intensity is a function of wavelength A, 
the spatial point r, the direction fi and the time. We want 
to restrict this to a simpler situation, where the intensity 
can be calculated from a plane parallel model atmosphere 
at each point. Furthermore we assume, that the surface 
structure is equivalent to a static model atmosphere down 
to optical depths larger than 1 for each time step. In this 
case the specific intensity is only a function of the wave- 
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<) 



length, the local flux F and the cosine between fl and the 
vector normal to the surface n 

h = h(F,fi), 

with /i = n • f2. In this approximation the dependence 
of the specific intensity on the surface element specified 
by the angles 9 and <fi and the time is not explicit but 
arises from the the dependence of the local flux on these 
quantities. 

The surface coordinates (6, 4>) can be taken as iden- 
tical to the coordinates in the spherical harmonic of the 
pulsation mode. For the integrated flux it is convenient 
to introduce a coordinate system with the axis oriented 
towards the observer. The system ip) is defined by 

d = for /i = 1. 

In this case the integrated flux can be written as 

p2n pi 

Fx(n,t)= dip dcostf cos& I x (F(<&,<p,t),iJ,). (11) 
Jo Jo 

Assuming that the flux J-\ is strictly a periodic function 
with the circular frequency u> we can expand the flux in a 
Fourier series 



J-\(fl, t) = + [An cos ncot + B n sinnujt] 

71=1 

with 



A n — 



and 



t +2-K/u> 



dt!F\(£l, t) cos ncot 



to 



B n = — / dtJ-\(Q,t) sinnujt. 

* J to 



(12) 



(13) 



Aq/2 is the mean surface- integrated flux. The amplitude 
for each harmonic n is defined by 



6T_ _ + Bl 

T n A 

and the normalized amplitude spectrum by 



Qn = 



-TnW A (X ) ^Al{\) + Bl{\) 



f n (Ao) A {X) ^42 (Ao)+i? 2 (Ao) 



(14) 



(15) 



Here Ao is an arbitrary reference wavelength, often taken 
in the visual at 5500 A corresponding to the V magnitude. 



5.1. Relation of I\ to the local flux F 

For the evaluation of Eq. (^lj) it is necessary to connect 
the specific intensity with the total flux of a local column. 
Using an expansion to first order we assume for this rela- 
tion 



(jx)F(t) (16) 



Expanding the local flux in a Fourier series in the same 
manner as for the surface integrated flux 



F(t) 



«o 



[a n cos nut + b n sin nut] 



the amplitude for the harmonic n is 



SF _ 2^fc 



bl 



F 



ao 



(17) 



Inserting (|16j) into (|ll| ) one obtains after some manip- 
ulation (Wu 1998 ) for n > the very simple result 



A„ 



B„ 



2 \ l\ , 
dip / dcosv cosv 

n Jo 



dF 
Oh 



2tt pi 

dip d cos •& cos ■& „ 
o Jo dF 



((!(■&, ip))a r , 



(fj,(0,ip))br, 



(18) 



(19) 



The integration for Ao leads to two additional terms 
and reduces to the simple form above only if the local 
time-averaged flux ao/2 is equal to the equilibrium flux 
Fo for the local column. 

Using approximations similar to those of Robinson et 
al. (1982) we could simplify the relation for the intensity 
as 



Ix = h x (n)F 



(20) 



with a limb darkening function hx- The major simplifi- 
cation is that in this form the limb darkening does not 
depend on the total flux and the intensity derivative in 
Eq.([l6]) is simply given by h x- In this case the integral for 
the global coefficient Ao always reduces to the simple form 
of Eq.([l8]). However, we will show below that this approx- 
imation is not adequate, since the the amplitude spectra 
depend very sensitively on the limb darkening and in the 
derivative of Eq. ( fj()|) 



dlx 
dF 



hx(v,Fo) + 



(21) 



Fo 



the second term is in general not negligible. 



5.2. The linear assumption 

The linear assumption is to specify the local flux by 



F 



SF 1 
F i 



17 



(22) 



If we want to determine the global flux amplitudes di- 
rectly from the local amplitudes without the intermediate 
step of Fourier coefficients we need to make further as- 
sumptions. If the coefficient Oo and the ratio a n /b n are 
constant over the stellar surface, it is possible to change 
the order of the square root and the integration over the 
visible disk in Eqs.( pT|Jl^ , [l9| ). After some algebraic manip- 
ulations we get 



f 



ST 

7, 



^ dp £ dcos d cos d ^\ Fo (M^ V)) SJ ' 



So" dl ? $0 dcosd C0S§ ^f\fI^^^ 



(23) 
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Fig. 9. Comparison of the amplitude spectra for an I = 
2 mode, calculated with a limb darkening independent of 
T ff (solid line) and with the linear expansion of the inten- 
sity (dashed line) . The spectra for the I — 1 mode are very 
similar and represented by the dotted line in the plot. 

It is common to transform the spherical harmonic 
Y l m (6,(j)) to a system of spherical harmonics in the co- 
ordinate system ip) by 

YT(M)= E R l m >mYi m X$M 

-Km'<l 



As was shown in several papers (e.g. Robinson et al. 1982) 
all terms with m' ^ do not contribute to a flux variation 
in the integrated light, because they are equivalent to a 
rotating pattern with respect to an axis directed to the 
observer. For the term with ml = the assumption of 
a constant phase over the stellar surface is valid and the 
Fourier coefficient a is assumed to be constant over the 
surface by Eq.(^2|) implicitly. Consequently, we can use 
Eq. ( p3|) to calculate the amplitude of the integrated flux 
with the local amplitude 



S_F' 
F i 



R, 



0m 



5F_" 



JViP,(cOB1?). 



(24) 



The quantity Ni symbolizes the normalization of the 
spherical harmonic YJ° and P; is the Legendre polyno- 
mial. The only inclination dependent term in Eq. (p4|) is 
the coefficient R l 0m . The amplitude {SF/F)^ X and R l 0m 
are constant over the whole surface and wavelength inde- 
pendent. Consequently, the normalized spectra Q are in- 
dependent of amplitude and inclination. Fig.^ shows one 
example of the results for a flux independent limb dark- 
ening and for the more general form of the linearized in- 
tensity for the model parameters we use in this paper, but 
with the linear assumption of Eq.(p^). Both results are 
identical for an I = 1 mode, but differ for larger I. Even 
for the I = 2 modes, there are significant differences. 

To test the accuracy of the linear approximation for the 
intensity (Eq.(|l6|)) we have calculated the Fourier coeffi- 
cients directly from the integrated light curve by means 
of Eq.(|l"l"|), ( |l2|) and ([l3| ) and have compared the result 
with the results from the linearized intensity. We find a 
very good correspondence for all wavelengths for modes 
with I < 3 and amplitudes less than 30%. This range is 



much larger than the range, where the total flux variation 
can be assumed to be linear with T e ff . For even larger am- 
plitudes the amplitude spectrum calculated directly from 
Eq.(^) becomes inclination and amplitude dependent, be- 
cause the function h\ varies over the surface for different 
flux amplitudes. Such large amplitudes are not observed, 
however. 

5.3. The Goldreich and Wu theory 

The Goldreich and Wu theory assumes - - following 
Brickhill — a pressure variation with a spatial and time- 
dependence described by the spherical harmonic and pe- 
riod of the pulsation mode. They predict the appearance 
of non-sinusoidal flux variations at the surface 

SF a SF . , . SF . .„ 

— = — + — sin(w*-ai) + — sm(2wt-a 2 ) + -(25) 
r I r i r 2 

An important result of their calculation is that the am- 
plitude of the fundamental is proportional to the pressure 
amplitude and the spatial distribution therefore still given 
by the same spherical harmonic. The first overtone varies 
quadratically with the pressure, which leads to the ap- 
pearance of a limited number of spherical harmonics. The 
phases a n are independent of the amplitude of the varia- 
tion. 

The result for the fundamental is equivalent to the 
linear assumption Eq.(p2|) and as a consequence leads to 
the same results for the normalized spectra: they depend 
only on the mode number I. 

The amplitude of the first overtone has a different be- 
havior than the spherical harmonic of the mode over the 
stellar surface. It can be written as 



SF _ SF ma * 

~F~2 ~ ~F~2 



[NiPi(cos6)e la2 



SF_ m: 

~F~2 



V{6,<t>) (26) 



The angle dependent function V can be expanded in a 
series of spherical harmonics 



(27) 



i'=0 i 



It is obvious that in this case V is a polynomial of the 
order 21 and consequently all coefficients c™ vanish for 
V > 21, but for this discussion we only need, that the sum 
has in general more than one term. Analogous to the linear 
assumption we can transform each term of Eq. ( p7[ ) to the 
coordinate system of the observer 



r>l ym 
J - L m"m' * V 



(0,<p) (28) 



1=0 ' 



-—v 



-—V 



The relation Eq. ( p8| ) can be introduced in Eq. (|2^) and we 
get an expression for (SF/F)2 in Eq.(p3|). Again, only the 
terms with m" = lead to a non vanishing amplitude 
in the integrated light. Inclination dependent quantities 
are only the coefficients R l 0m i, but they are different for 
all I'. Consequently, the inclination dependence does not 
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cancel in the normalized spectrum Q2 and the spectrum 
becomes inclination dependent. In contrast to this, the 
amplitude (5F/F)™ & * cancels out, and the normalized am- 
plitude spectrum remains independent of the amplitude. 

The 2 nd overtone will not be discussed in detail. The 
local flux amplitude (SF/F)^ cannot be represented by a 
single (<5F/F) max and consequently the normalized ampli- 
tude spectra become in general also amplitude dependent. 



6. Amplitude spectra from numerical calculations 

To determine the normalized amplitude spectra from our 
numerical results, it is possible to use the integration of 
the local Fourier coefficients a n and b n . However, because 
the phase is not constant (see Fig. ||) over the surface 
and do is a function of the coordinates (9, (f>) (see Fig. ||) 
the expression Eq. (^3|) is completely useless even for m = 
modes and we have to calculate the coefficients A n and 
B n separately. This means that there is no economical 
advantage to calculate the Fourier series locally and then 
to integrate the coefficients over the visible disk. 

We therefore decided to calculate the flux T\ for each 
time step completely numerically by using Eq. (pT|). This 
method also avoids the disadvantage of having to use a 
linearization like Eq. (|l^) to reach expressions for the 
Fourier coefficients of the integrated flux T\. From the 
time-dependent flux we can calculate the Fourier coeffi- 
cient with the aid of Eq. (|l2] ) and Eq. (|l3|) . The numerical 
results presented in the following section are based on such 
integrations. 

On the other hand, in order to gain more insight into 
the qualitative behavior of the spectra for different am- 
plitudes and inclinations, we will use simple fits for the 
numerical results of the local variation, which allow then 
analytic integrations for the total flux. Both approaches 
are complementary: with the completely numerical simu- 
lation we avoid simplifications, which may not be correct, 
but we can calculate the result only for very few combi- 
nations of parameters. The semi-analytical study, on the 
other hand, shows the functional dependence of the results 
on quantum numbers, inclination, etc. 



6.1. General effects on the spectra 

The phase dependence of all local Fourier coefficients and 
the nonlinear dependence on the pressure variation leads 
to deviations from the assumptions made traditionally for 
the flux variation over the stellar surface. The intention 
of this subsection is to show, how this deviation affects 
the synthetic spectra for different inclinations, pulsation 
amplitudes and mode numbers I. 

To find an expression for the error made when using 
the linear assumption, we assume a flux independent limb 
darkening function and the integrated Fourier coefficients 
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Fig. 10. The amplitude dependence of the correction 
term Q™ rr . The factor T is displayed for our model pa- 
rameters. 

given by Eqs. and ©. We define a correction term 
for the fundamental as 

¥ W ¥ (A) 

Qicorr J- n v ' J~ n v ' 



(29) 



linear 



In appendix A we demonstrate that in a first order ap- 
proximation Q™ rr can be separated into three factors 



QT 



A (X )I (X) / 



A (X) 



5P 

17, 



A(i,i). 



(30) 



The first term is practically independent of wavelength 
(To is defined in Eq. |35]in the Appendix). The second term 
r is derived from the expansion coefficients of the local flux 
as a function of the pressure variation. For finite ampli- 
tudes we interpret these coefficients as fit parameters to 
describe the local flux as a polynomial of 2 nd order up to 
the maximum pressure variation, that occurs at the stellar 
surface. This is plotted in Fig. (|l0|). T is not a function 
of the inclination or the mode numbers, but depends on 
the stellar parameter as T e ff and the pulsation period. For 
our model parameters we can divide the synthetic spec- 
tra in three different pressure regimes: Small amplitudes 
up to 5P/P m ax ~ 5%, where the function T is positive; 
intermediate amplitudes 5% < SP/P max < 7%, where T 
is approximately 0, and large amplitudes SP/P max > 7% 
with a negative T. In the following 3 subsections these 
cases are discussed separately for the numerical results. 

The third term A is the inclination dependent term. 
Besides this it depends on the mode numbers and wave- 
length, but not on the amplitude. The results for the factor 
A for I < 3 are plotted in Fig. (pT|). This factor is for most 
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Fig. 11. The factor A is shows the dependence of the 
correction factor between the linear assumption and the 
numerical simulation on inclination. The 3 curves show 
the factor for 3 different mode numbers (I = 1,2,3) and 
the Edington limb darkening law. Large deviations occur 
near an inclination of 55 degrees. 

inclinations very small for I = 1 and moderate for Z = 2, if 
the inclination is clearly different from the latitude of the 
node lines. For I — 3 the factor A is much larger and leads 
to a large global deviation of the numerical results from 
those using the linear assumptions for any inclination and 
amplitude. 

6.1.1. Small amplitudes 

In the previous section we have discussed the effects of the 
nonlinear behavior of the light curves on the amplitude 
spectra of the fundamental with some approximations. To 
find the correct spectra we now calculate the Fourier coeffi- 
cients directly from the integrated light curves in a strictly 
numerical way. 

As an illustrating example we take the I = 2 mode with 
a maximum pressure amplitude of 2%. For this amplitude 
the light curves are sinusoidal in a good approximation, 
so we can restrict the discussion to the fundamental of the 
mode. 

From Fig. ( Jjj~T| ) we do not expect very large deviations 
from the linear result for most inclinations. To show the 
effect we take an inclination of 60° near the latitude of 
the node line, where we expect a moderate deviation. The 
numerical result is plotted in Fig. (|l2j), together with the 
result of the linear assumption, that is supported by the 
theory of Goldreich and Wu. All amplitudes in this and 
the following similar figures have been normalized to 1 
at 5500 A. Although the effect is not very large, it is 



Fig. 12. Amplitude spectra for small pressure amplitude. 
The solid line and the dashed lines are results based on 
the spherical harmonics I = 2 and I = 1 respectively. The 
diamonds connected with the dotted line arc the numerical 
results. 

important for mode identification, because of the gener- 
ally small differences between the calculations for different 
mode numbers. 

Deviations in the spectra already for small amplitudes 
are suggested by the analytical calculations in the previous 
subsection, but are not obvious in the plot of the absolute 
amplitude Fig. (|§|) , which grows linearly with the pressure 
amplitudes up to ~ 5% pressure amplitude. Consequently, 
the deviation has to be an effect of the phase shift for 
small amplitudes. This phase shift should be visible in the 
integrated light as well. Fig. ([l3]) gives the numerical pre- 
diction for the phase shift for the same mode and modes 
with / = 1. In contrast to the prediction of Goldreich and 
Wu, we expect a small phase shift of a few degrees over 
the spectral range for the I — 2 modes. 



6.1.2. Intermediate amplitudes 

The amplitude dependent factor T is very small in this 
region and we expect practically no differences to the pre- 
dictions of Goldreich and Wu. Fig. (14) shows the numer- 
ical result for the same mode and inclination as Fig. (|l2|), 
but for a pressure amplitude of 5%. The spectra for the 
numerical prediction and the Goldreich and Wu theory 
are very similar. 

In contrast to the small amplitude case, the light 
curves for 5% pressure amplitude are non-sinusoidal with 
large flux amplitudes. We can calculate the expected devi- 
ation from the predictions of Goldreich and Wu analogous 
to Eq. ( fiol) and find for this region a small deviation for 
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Fig. 13. Phase spectra for small pressure amplitude. The Fig. 15. Amplitude spectra for the first overtone and in- 
solid line (I — 1) and the dashed line (I = 2) gives the pre- termediate amplitude. The lines have the same meaning 
dieted phase shift for the wavelength range from 2000 Ato as in Fig (|l2|). 
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Fig. 16. Amplitude spectra for large amplitude. The lines 
Fig. 14. Amplitude spectra for intermediate amplitude, have the same meaning as in Fig @. 
The lines have the same meaning as in Fig (E2I). 



6.1.3. Large amplitudes 

the first overtone as well. Fig. dl5|) shows a comparison 

r , , . , .,, N — 1 , , , r For increasing maximum pressure amplitude the devia- 

ot the numerical result, with the predicted spectrum tor a . , r . 

, , . , , r ,1 a ,., , ,1 tion factor T becomes smaller and finally negative. For 

quadratic dependence of tnc flux amplitude on trie pres- 

,., j really large pressure amplitude we expect again a signifi- 

sure amplitude. , r jL3 , , , (3 1 to i — 1 

cant deviation from the Goldreich and Wu result. Fig (|16|) 



shows the numerical result for a 15% pressure amplitude. 
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As expected from the analytical result, the difference to 
the Goldreich and Wu result has the opposite sign and the 
spectrum of the I = 2 mode becomes very similar to that 
of an I = 1 mode. The numerical result for the / = 1 mode 
does not significantly change for an inclination of 60° (see 
Fig [ll]) and modes with I = 2 can easily be confused with 
1 = 1. 
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Fig. 17. Amplitude spectra for the first overtone and 
large amplitudes. The lines have the same meaning as in 
Fig ©. 



The first overtone is significant for large amplitudes 
as well. In contrast to the intermediate amplitudes, a dif- 
ference from the predicted result of Goldreich and Wu 
appears. For our example Fig. ( Jl7| ) shows the numerical 
result in comparison with the quadratic dependence of the 
flux amplitude. 

For I = 3 the deviation of the numerical results 
from the results of the linear assumption is much larger. 
Fig. ( |l8| ) shows one example for a large deviation. In con- 
trast to a rapidly increasing amplitude in the UV the am- 
plitude remains almost constant. In this regime the pre- 
diction of the variation from Eq. ( |30| ) can be only quali- 
tative. In general the amplitude increase towards the UV 
is caused by the decreasing "effective visible area" due to 
the larger limb darkening. This effect becomes significant, 
when the effective area is comparable in size to the spatial 
scale of the light variations on the surface of the star. For 
the large amplitudes the spatial structure of the I = 3 flux 
variation is more complex than the unperturbed spherical 
harmonic and the cancellation effects set in only farther 
in the UV. 
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Fig. 18. Amplitude spectra for I = 3 and large ampli- 
tudes. The lines have the same meaning as in Fig [l2[ 



7. Conclusions 

In the past the method of time-resolved spectroscopy has 
been based on the assumption — called "linear" through- 
out this paper — that the flux varies with the spherical 
harmonic of the pulsation mode over the surface of the 
ZZ Ceti stars. The numerical light curve simulation shows 
that this assumption is questionable in many situations. 
For small amplitudes with sinusoidal lightcurves, only the 
absolute amplitude varies like the spherical harmonic of 
the mode. The phase shift leads to normalized amplitude 
spectra, which depend on inclination and pulsation (pres- 
sure) amplitude for mode numbers I > 2. For large am- 
plitudes (> 10% in pressure) the absolute amplitude of 
the flux depends non-monotonously on the pressure am- 
plitude and the flux variation shows maxima at different 
latitudes than the spherical harmonic of the mode. The 
numerical results show relatively good correspondence to 
traditionally calculated amplitude spectra for intermedi- 
ate amplitudes (~ 5%) in pressure, with non-sinusoidal 
light curves. In this regime the phase becomes stationary 
and the absolute amplitude varies in good correspondence 
to the spherical harmonic of the mode; the numerical re- 
sults are close to the results of the perturbation analysis 
of Goldreich and Wu. 

Another result of our simulations is the existence of 
a maximum surface flux amplitude. This is a reaction of 
the convection zone on a predefined pressure variation in 
deeper layers, and should not be confused with the ampli- 
tude satur ation of the pulsation as studied by Goldreich & 
Wu (1992). The amplitude of the photospheric flux is re- 
duced by the inert reaction of the temperature structure 
in the convective layer to the changing input flux. The 
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decrease in the thermal heat flux is converted to kinetic 
energy of the pulsation as discussed above, and this con- 
vection region is therefore a driving region for the pulsa- 
tion. The extent of this conversion depends on the thermal 
time scale of the convection zone. 

In our model calculations for large amplitudes the 
time-averaged depth of the convection zone becomes much 
larger than in the corresponding static model. The thermal 
time scale then increases with the pulsation amplitude, 
leading to a reduction of the surface flux amplitude. This 
is the reason for the existence of a maximum amplitude for 
the surface flux. It also explains the return to sinusoidal 
variations for large pressure amplitudes: the convection 
zone becomes so thick that variations during one cycle do 
not decrease its depth enough for the appearance of higher 
overtones. 

While our calculations cannot determine the maximum 
pressure amplitude that can be reached in a pulsating star, 
one prediction is that observable flux amplitudes should 
not exceed 10%, for the parameters used in this paper. The 
maximum flux amplitudes support the dominance of small 
I in the observed light curves and lead to the prediction, 
that the dominant pulsation modes are I = 1, if / = 1 
modes and modes with larger I are unstable, for a large 
amplitude range, independent of the actual amplitudes of 
the modes. 

While nonlinear effects are well known for large ob- 
served amplitude pulsations, the present study reaches the 
surprising conclusion, that the amplitudes and phases can 
be described fairly accurately by the linear theory as de- 



veloped by Robinson et al. (1982) for the fundamental and 
Wu ( 1998| ) for higher harmonics. 

On the other hand, the observation of small ampli- 
tudes in the surface integrated light does not necessarily 
guarantee that these effects are absent. Deviations from 
the traditional interpretation could occur in spite of small 
observed amplitudes under the following conditions 

- the pressure amplitude is larger than about 12%, that 
is in the range were the flux amplitude decreases again 
with increasing pressure amplitude 

- the flux variations are intrinsically large on the surface 
but the inclination is close to the latitude of a node line 

- the mode number I is large 

In all these cases one of the correction factors describing 
the difference between our numerical results and the am- 
plitudes calculated with the linear assumptions becomes 
large, and mode identification from time-resolved spec- 
troscopy will only be possible with extended numerical 
simulations. 

All conclusions are obtained for a special set of stellar 
and pulsational parameters. The effects discussed may be 
more or less important for other ranges of log g and T e g. 
We have also not taken into account that in general more 
than one pulsation mode is present in a star. The modes 
will be influenced by the presence of other modes and 
by their properties even well below the convection zone. 
These effects are beyond the scope of the present analysis. 



We do not, however, expect that the main conclusion will 
change: the mechanism that leads to non-sinusoidal light 
curves for large flux variations can be important for a 
correct identification of the pulsation modes using time- 
resolved spectroscopy, for any observed amplitude of the 
flux. 
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Appendix A: Difference between the linear spectra 
and the numerical results 

To find a semi-analytic expression for the difference be- 
tween the synthetic amplitude spectra from the linear the- 
ory and our numerical calculation, we approximate the 
local Fourier coefficients a\ and b\ as a quadratic func- 
tion of the local pressure amplitude. Assuming a pressure 
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variation with the spherical harmonic of the mode on the 
stellar surface we can write for the modes with I = 



ai « cti 



6P_ 

P max 



cv 2 



5P 

P max 



(Pi? 



(31) 



The quatities a.\, ct2 are given by the fit to the numerical 
results and thus depend on the maximum pressure ampli- 
tude on the surface. 

In a similar way as for the spherical harmonics we can 
also expend (P;) 2 with a finit number of P;/. In this first 
order approximation we only need the term with I' = 
of this expansion (the next higher term belongs to I' = 2 
and is strongly reduced by the cancelation effect). We can 
write 



(PO 2 



1 



(32) 
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Now it is possible to perform the rotation to the coordi- 
nate system of the observer in the space of the spherical 
harmonic. The result is a mixing of all possible quantum 
numbers m! to the same I. As in the linear case we only 
have to consider the contribution of m' = to the visible 
pulsation amplitude and get 

a\ w aiP(cosz) P(cos^) + a 2 ~^j~~£ > 

or after the integration over the visible disc 

1 
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To find the Fourier amplitude we insert this expression 
and the anlogous equation for B\ in (|l4|) and get for small 
oc2,02 (Pi and (32 are the fit parameters for B\) 
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or by expanding the sqare root 
ST 2 
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The first term of this expression is proportional to the 
result of the linear theory with the assumption of a flux 
variation with a spherical harmonic. The second term is 
the first order nonlinear correction term for the numerical 
result. For this first order approximation it is sufficient 
to identify the mean flux An/ 2 with the flux of the non- 
pulsating star Jo. We get for the correction term of the 
amplitude 

Fx * ' (2l + l)^+~M 



(38) 



This expression is a function of T e g of the equilibrium 
model, the pulsation mode and the maximum pulsation 
amplitude of the model; it does not depend on the wave- 
length nor on the inclination. The linear term, however, 
which is normally dominant, is a very sensitive function 
of the inclination. 

Following the standard strategy to find an expression, 
which eliminates the inclination dependence of the linear 
term, we normalize Eq. ( pTj ) to the amplitude at a refer- 
ence wavelength Ao (e.g. = 5500 A) and obtain for the 
normalized spectrum 

Qi w Qi in + Qi orr 

with the linear normalized spectrum Q\ ln and 



r .WW r (f)i W 

The inclination dependent term A(i,l) is given by 

w. .v = __}- 1 

(h '- 2i + lP,(coBt)Zi(Ao) 
and the amplitude dependent term by 
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